An Efficient Algorithm for Density Functional Theory Simulation of Large Quantum Dot Systems 
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Kohn-Sham spin-density functional theory provides an efficient and accurate model to study electron-electron 
interaction effects in quantum dots, but its application to large systems is a challenge. An efficient algorithm for 
the density-functional theory simulation of quantum dots is developed, which includes the particle-in-the-box 
representation of the Kohn-Sham orbitals, an efficient conjugate gradient method to directly minimize the total 
energy, a Fourier convolution approach for the calculation of the Hartree potential, and a simplified multi-grid 
technique to accelerate the convergence. The new algorithm is tested in a 2D model system. Using this new 
algorithm, numerical studies of large quantum dots with several hundred electrons become computationally 
affordable. 
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PACS numbers: 

I. INTRODUCTION 

Quantum dots (QD) are one kind of nano-device in which 
the motion of electrons is quantized in all three dimen- 
sions through the lateral confinement of a high-mobility 
modulation-doped two-dimensional electron gas in a semi- 
conductor hetero-structure. '-^-^ Both quantum interference 
and electron-electron interactions play important roles, and 
their interplay underlies many properties that are fascinat- 
ing in terms of both fundamental mesoscopic physics and fu- 
ture technical applications. Various theoretical models have 
been developed to explain experimental discoveries and to 
predict new properties. The Kohn-Sham (KS) density func- 
tional theory (DFT) method"'-^ '' provides an accurate numeri- 
cal model to study electron-electron interaction effects in QD 
systemsii2i2ii2iiiii2iyiiiii^ii^iiLi^ In spite of its comparatively 
low computational cost, previous DFT calculations are limited 
to systems in which the electron number is generally less than 
a few tens. In many experimental cases, however, the elec- 
tron numbers involved are more than several hundred. The 
main theoretical approaches in the regime of large number 
of electrons are statistical methods^i^iiSiSS which are usually 
based on some general assumptions whose validity, in many 
cases, is yet to be justified. It is therefore desirable to explore 
the large dot regime directly by using numerically more ac- 
curate models such as DFT. This imposes a demanding com- 
putational task because to obtain meaningful statistics, many 
calculations with several hundred electrons need to be done. 
It is therefore compelling to develop more efficient numeri- 
cal techniques for DFT simulation of QD systems; this is the 
main object of the current study. 

Two key issues are involved in the numerical implemen- 
tation of the KS-DFT methodi^ii^^ (1) the numerical repre- 
sentation of wave functions and the KS Hamiltonian, and (2) 
the solution of the numerical KS equation. While local ba- 
sis sets (mainly Gaussian-type orbitals) dominate in conven- 
tional quantum chemistry of molecular systems, the plane 
wave (PW) basis, ^''^^ combined with the pseudo-potential 
method, is widely used in the ab initio electronic calculations 
of various material systems, in which the fast-Fourier trans- 



form (FFT) method can be used to take full advantage of the 
periodicity of crystal structures. In principle, the PW method 
is valid only for periodic systems, but aperiodic systems can 
be treated by introducing the super-cell technique.^' In recent 
years, several groups have been advocating the use of basis- 
free real space methods for electronic structure calculations of 
finite systems, in which the wave functions are represented in 
real space, and the kinetic energy operator is discretized by a 
high-order finite difference (FD) method. With a given rep- 
resentation, there are various methods to solve the resultant 
numerical KS equation. Roughly they fall into two different 
types: methods that minimize the total energy directly, and 
those that solve the KS equation in a self-consistent way. In 
addition, for DFT simulations of finite systems, another im- 
portant issue is the calculation of the Hartree potential. 

Aiming at modelling QD systems efficiently, we have de- 
veloped new techniques in both the numerical representation 
of KS orbitals and the solution of the KS equation. We note 
that in QD systems, the wave functions vanish at the bound- 
aries, and in some cases even the hard-wall boundary condi- 
tion is used. For a function in a rectangular box with zero 
boundary values, the most natural basis set is the particle-in- 
the-box (PiB) basis set. The kinetic energy operator is di- 
agonal in the PiB space and the transformation between real 
and PiB space can be efficiently performed by the fast sine- 
transform (FST) method, which is a variant of the FFT. For 
the solution of the KS equation, we modified Teter, Payne and 
Allan (TPA)'s band-by-band conjugate gradient method^i to 
get a more efficient direct minimization approach. 

The outline of the paper is as follows. In the next sec- 
tion, after a simple description of the KS-DFT method, we 
present the main components of our algorithm: (1) the PiB 
representation of the KS equation; (2) a modified band-by- 
band conjugate gradient method for the direct minimization 
of the KS total energy; (3) a Fourier convolution method for 
the calculation of the Hartree potential that was developed by 
Martyna and Tuckerman^^ in the PW pseudo-potential calcu- 
lations; and (4) a simplified multi-grid technique to accelerate 
the convergence. In Section III, the new algorithm is tested in 
a 2D model QD system with electron number N = 100. Sec- 
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tion IV summarizes the main results and concludes the paper. 

II. METHOD 
A. Kohn-Sham spin-density-functional theory 



approximation (GGA) are available, it has been shown that the 
GGA results are close to those from LSDA calculations in QD 
systemsiii In terms of the implementation, the calculation of 
Vxc is trivial when the spin densities are in real space as is the 
case in our algorithm, but the calculation of Vh requires more 
efforts as will be shown later. 



Considering the important role played by electron 
QD systems, we take the effect of spin polarization 
itly into account in the framework of Kohn-Sham spin- 
functional theory (KS-SDFT).'*'' In KS-SDFT, the 
state energy of an interacting system with electron 
N and the total spin S in the local external potential 
is written as a functional of spin densities n'^ with a 
denoting spin-up and spin-down, respectively. 



spm m 
explic- 
-density 
ground 
number 
Kxt(r) 
= a,l3 



E[n",n'^] =Ts [n",n^] + J 



n{r)Vcxtir)dr 



n(r)n(r') 



(1) 



(Effective atomic units are used through the paper: for GaAs- 
AlGaAs QDs, the values are 10.08 meV for energy and 10.95 
nm for length.) is the kinetic energy of the KS 

non-interacting reference system which has the same ground 
state spin density as the interacting one, and E^c [n" , n^] is 
the exchange-correlation energy functional The spin densities 
n" satisfy the constraint / n'^{r)dr = N'^ with N" = [N + 
2S')/2and iV^* = (iV-2S')/2. 

Assuming that the ground state of the non-interacting refer- 
ence system is non-degenerate, the non-interacting kinetic en- 
ergy is given by [n",n^'] = Y.t,a i^i 
the ground state spin density is uniquely expressed as 



(2) 



Here tpf are the lowest single-particle orbitals which are ob- 
tained from 



H^gVnr)=e7Vr(r), 
with the KS Hamiltonian Hks defined as 

H^s = -^V' + ^-t(r) + VH[n;v] + V:,[n-.u 



0.. 



(3) 



(4) 



yff[7i;r] and V^^[n" , n'^ ; r] ai'e the Hartree and exchange- 
correlation potentials, respectively. 



VH[n;r] 
K''Jn",n'3;r] 



n(r') 



d^r' 

r — r'l 



(5) 
(6) 



We have used the local spin-density approximation 
(LSDA)ii& for E^C7 which is widely used for the modelling 
of material systems. Although more accurate exchange- 
correlation functional forms such as the generalized-gradient 



B. PiB Representation 

To simplify the notation, we will take 1-D systems as an 
example, the generalization to higher dimensional cases being 
straightforward. Any regular function f{x) that is localized in 
the finite region < x < L with zero boundary values can be 
expanded as 



2 . nnx 



(7) 



and the expansion coefficients C„ are 



£ f{x)su,^dx. (8) 



Integrating Eq. (|8} numerically on a set of equally spaced 
discrete points {xj = jA: 
zoidal formula2& leads to 



discrete points {xj = jAx = j^} using the extended trape- 



a = ^/|ax''£' /, sin ^ = (9) 



L 

with/j = fixj) and 



i=i 



. njn 
3in 



= FST{/,} (10) 



where FST{/j} denotes the fast sine-transform of the data 
{/.}■ 

One of the key ingredients in our algorithm is the action of 
the single-particle Hamiltonian operator on wave functions. 



Uf(x)^Tf{x)+Yf{x) 



(11) 



the efficiency of which is critical for the performance of the 
whole algorithm. Whereas the potential energy operator is di- 
agonal in real space, the kinetic energy operator is diagonal in 
PiB space. Wave functions can be transformed between the 
two spaces efficiently by the fast sine transform. In our algo- 
rithm, wave functions are in discrete real space {fj}, and the 
application of the potential energy operator is therefore triv- 
ial, V{/j} ~ {Vjlj}, where Vj is the value of the potential 
at the point xj. The kinetic energy operator is applied to wave 
functions in PiB space. 



T/(a;) 



2m 



V2/(a;) = — Vi^n^sin— , (12) 
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with kn = mr/L, which in the discrete form becomes 



where 



2m 



}• 



(13) 



The PiB representation formulated above is closely related 
to the PW method. In fact, for finite systems with soft-wall 
boundaries, the two representations are numerically equiva- 
lent. Mathematically, however, they are different: The PiB 
basis set is real and the zero boundary condition is imposed by 
the basis set itself, but in the PW case the basis functions are 
complex and the wave functions are forced to be zero at the 
boundaries by the external potential of systems under study. 
The PW method will fail in the case of the hard-wall bound- 
ary problem, which is quite common in the studies of quantum 
dots, but the PiB method is still valid. 



C. Direct-minimization conjugate gradient metliod for the 
Kohn-Sliam equation 

In direct minimization approaches, the KS total energy 
functional is minimized directly over orbital wave functions 
under orthonormal constraints. We have made several impor- 
tant modifications to TPA's band-by-band conjugate-gradient 
scheme to obtain higher efficiency. Here we give an outline 
of the algorithm, and emphasize the modifications we have 
made. 

The basic idea of a band-by-band scheme is to minimize 
the total energy over one band (or orbital) at a time, which, 
compared to other conjugate-gradient schemes, ^^'^^^ has the 
following advantages: (1) much lower requirement for stor- 
age space, (2) simpler implementation, and (3) particularly 
in our algorithm, an efficient approximate line-minimization 
scheme, as will be shown. 

First the steepest descent (SD) vector for the i-th orbital at 
the m-th iteration is calculated from 




with A™ = (V'l"! Hks \4'T') ('^o simplify the notation, we use 
Dirac's state vector notation and drop the spin index in the 
following formulation). In TPA's algorithm, the SD vector is 
preconditioned before it is used to build the conjugate vec- 
tor. In our calculations, however, it was found that although in 
many cases preconditioning does accelerate the convergence, 
its effect is not always positive. On the other hand, the compu- 
tational overhead in the preconditioning step, which involves 
another orthogonalization process as well as the action of the 
preconditioning operator on the SD vector, can be expensive 
for large systems. As shown later, by using a simplified multi- 
grid technique, we can achieve fast convergence even without 
preconditioning of the SD vectors. 

The conjugate vector |<p™) is then constructed as a linear 
combination of the SD vector and the previous conjugate 
vector. 



(16) 



with 'yI = 0. The conjugate vector is further orthogonalized 
to the present band | V"™) and normalized (N is denoted as the 
normalization operator). 



Na-iv-rx^rDi^r 



(17) 



The new wave function for the i-th orbital li/"!"^^) formed 
from the linear combination 



|(^r>sin0min 



(18) 



which is guaranteed to remain normalized and orthogonal to 
all other orbitals. ^min is obtained by minimizing the to- 
tal energy as a function of 9 with \^i{9)) = \ip"^) cos9 + 
\(p'"^) sin 9. In TPA's algorithm,^''^"^ 6'niin is determined by the 
following approximate scheme: The total energy as a function 
of 9 is approximated by E{9) w Eavg+Ai cos 29+Bi sin 29; 
the three unknowns, Eavg, A\ and B\, are determined accord- 
ing to three pieces of information: E{9 = 0), ^^^|e=o and 
E{9 = 7r/300). 

Here we propose a more efficient approximate scheme for 
the determination of ^min- The derivative of E{9) with respect 
to 9 can be obtained from 

^^^^^ =2{ip'r\HKs{9)\^Pr)cos29 



89 



((V^riHKsWIV-r) - ((/:'r|HKs(^)l¥'r»sin2^ 

(19) 



Assuming Hks(^) ~ Hks(O), from ^^p- = we get 

1 _i B 



tan 
2 A 



(20) 



Hks) W) (14) with 



A = - (v^ri hks(o) \^t) + (^n hks(o) wn (21) 



and 



i3 = 2(v.nHKs(0)|C 



(22) 



(15) 



The underlying approximation in our scheme is similar to that 
of TPA's, but our scheme is much more efficient in terms of 
computational effort: In TPA's scheme, at each band iteration 
the total energy must be calculated twice, i.e. E{9 = 0) and 
E{9 = 7r/300); in our scheme, the most time-consuming step 
is the action of Hks on \(p'"''), which is much faster than the 
calculation of the total energy. Considering further the fact 
that the total number of band iterations can be very large, we 
note that an efficient line-minimization scheme such as ours 
is crucial to reduce the computational effort. 

In TPA's algorithm, after the wave-function is updated ac- 
cording to Eq. ( I18> . the Kohn-Sham Hamiltonian is updated 
immediately, which involves the reconstruction of Vh (r) and 
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14c (r) according to the new density. This is actually quite ex- 
pensive for large systems. On the other hand, we expect that 
the KS Hamiltonian will not experience large changes inside 
the iterations of a single orbital. So in our algorithm, we up- 
date Hks after every A^updatc band iterations, and the optimal 
value of A'^updatc will be explored in the next section. 

In each orbital, the procedure described above is repeated 
-^band times; the iterations are then started on the next or- 
bital. After the wave functions of all orbitals are updated in 
this way, the total energy is calculated and is compared to that 
of the previous cycle to determine if the final convergence is 
achieved. The main parameters in the algorithm are A'^band 
and iVupdatc- Their effects on the performance of the algo- 
rithm will be tested in detail in the next section. 



D. Calculation of Vh 

For finite systems, the simplest and perhaps most ineffi- 
cient way to calculate Vh is by direct numerical integration, 
which is feasible only for small systems. Another widely- 
used approach is to solve the Poisson equation equivalent to 
Eq. (|5}- Though the Poisson equation itself can be solved 
with great efficiency, the calculation of boundary values can 
be quite expensive even by using efficient multipole expan- 
sion techniques. Additionally, the Poisson solver approach is 
valid only for 3D systems; in the case of 2D systems, there 
is no Poisson equation equivalent to Eq. (|5}- In recent years, 
several schemes have been proposed to extend the conven- 
tional Fourier convolution method to finite sy stems. ^.^i22i^2ii 
In particular we have incorporated Martyna and Tuckerman's 
method^^ into our algorithm. Considering that the Martyna- 
Tuckerman method was developed mainly for the modelling 
of molecular and material systems within the plane-wave 
pseudo-potential framework, we will formulate the approach 
here with some detail. 

The calculation of the Hartree potential is straightforward 
for periodic systems, but this is not the case for finite aperi- 
odic systems. The potential Vf/(r) has the form of the convo- 
lution between the density and the Coulomb interaction ker- 
nel, Vc{r) = 1/r, which has the following simple relation in 
the Fourier spacep^ 



yf^(k) =n(k)«,(k), 



(23) 



where /(k) refers to the Fourier transform of /(r). Eq. j23t is 
useful only when we have the analytical form of n(k) and can 
perform the inverse Fourier transform of Vh (k) analytically, 
which is not true for most cases where the density is usually 
represented in discrete real space. 

When applying the Fourier method to discrete finite sys- 
tems, periodic boundary conditions are always assumed. It 
has long been known that the unphysical interactions between 
neighboring super-cells can be avoided by calculating Vh in a 
doubly extended grid.^- In particular for 2D systems as illus- 
trated in Fig. n the original x Ly grid (il) is extended to 
2Lx X 2Ly (il2L)- The density in the extended grid is defined 





n^T (r)=n(r) 


O' 

■ # 


naL(r)=0 : 



Q 



2L 



FIG. 1: Illustration of the Fourier convolution method for the cal- 
culation of the Hartree potential in finite systems. The original grid 
where the density and potential are defined is extended into doubled 
space. The density in the extended region is set to zero. Imposing 
the periodic boundary condition for this extended grid, the unphysi- 
cal interaction between neighboring super-cells can be avoided. 



as 



n2L(r) 



_ / n(r) 




if r e 
otherwise 



(24) 



Imposing periodic boundary conditions to both the density 
and the Coulomb interaction kernel in the extended grid, the 
potential can be calculated according to the convolution theo- 
rem. 



yjf(r) =^n2L(k)i;e(k)e 



ikr 



(25) 



where n2L(k) and iJc(k) are respectively finite Fourier inte- 
grals of the density and the Coulomb interaction kernel in the 
extended grid. 



n2L(k) = 7^ / n2L{r)e '^''dr, 
zJc(k) = / i;,(r)e-*-^dr 



(26) 
(27) 



where VI2L is used to denote both the extend grid and its vol- 
ume (or area in the 2D case). 

While rT2L(k) can be easily obtained from the discrete 
Fourier transform of its real space values by EFT, the calcula- 
tion of TJc(k) is much more involved because of the singular- 
ity of the Coulomb interaction kernel in real space. The key to 
Martyna-Tuckerman's approach is to decompose the Coulomb 
interaction kernel into long- and short-range parts. 



(r)+w(^''°'-*)(r) (28) 



where erf {x) and erfc(a;) are the eiTor function and its com- 
plement, respectively, and a is the parameter that controls the 
effective cut-off range. The finite Fourier integral of the short 
range part can be well approximated by its infinite Fourier 
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transform, 

Jwhole space 

(29 

which is analytically known in both 2D and 3D cases. The 
finite Fourier integral of the long-range interaction can be di- 
rectly obtained from the discrete Fourier transform of its real 
space values. In the practical implementation, Vc needs to 
be calculated only once at the beginning. The calculation of 
Vff (r) involves only two FFTs (one forward and one back- 
ward), which makes this approach much more efficient than 
methods based on a Poisson solver. 




32 48 64 80 

Number of Grid Points in Each Direction 



FIG. 2: Convergence of the total energy with respect to the number 
of grid points using 5-Point FD, 13-Point FD and PiB representation. 



E. A simplified multi-grid technique 

The multi-grid method is an efficient technique to ac- 
celerate the convergence in various real space relaxation 
approaches. The basic idea is that low-frequency errors are 
easier to eliminate in a coarse grid than in a fine grid. Lee et 
al.^^ proposed a simple one-way multi-grid technique in finite- 
difference real space KS calculations. Wave functions being 
represented in real-space in our algorithm, a similar simplified 
multi-grid (SMG) method can be implemented in a straight- 
forward way: the KS energy functional is first minimized in a 
coarse grid; the converged wave functions, after interpolation 
and re-orthogonalization, are taken as the initial guess for the 
minimization on the fine grid. With these well-preconditioned 
initial wave functions, the convergence in the fine grid can be 
easily attained. 



III. NUMERICAL TESTS 

In most experimental QD systems, the excitation in the ver- 
tical direction can be neglected, and a 2D model is a good 
approximation. In this paper we will report test results only in 
a 2D system; the extension to 3D systems is straightforward, 
and the high efficiency of our algorithm makes it possible to 
explore large dots even in the 3D case. We test the perfor- 
mance of the new algorithm in a coupled quartic oscillator 
potential system. 



K=xt(r) 



+ by — 2Xx y +"f{x''y — xy )r 



(30) 

with a — 10^*, b — tt/4, A = 0.6 and 7 = O.L Considering 
that the convergence behavior in KS calculations is generally 
related to the effective interaction strength, which is charac- 
terized by rg = {nn)^^^^ /qb in the 2D case where n is the 
average density and as is the Bohr radius, we have chosen 
the parameters in I^xt so that the estimated is about L5, 
which is close to experimental values. The calculations are 
done in a grid of size = Ly = 50 and the number of grid 



points is Nx = Ny = 64. All the following numerical results 
come from calculations with electron number N = 100 and 
spin S — 0. For the exchange-correlation energy, E^c, we 
use Tanatar and Ceperley's parameterized form of the LSDA 
functional.^'* The convergence criterion is set as e = 10^^, 
which corresponds to about 10^^ meV in GaAs-AlGaAS QD 
systems. 

Considering that the FD method has been widely used in 
the numerical modelling of QD systems, ^-^Jili we first make 
a comparison between FD and PiB. In the FD representation, 
the second order derivative in the kinetic energy operator is 
locally discretized in real space 



1 

ft2 



J2 Q/, +0(/i2™) (31) 



j = -m 



where h is the discretization step and the coefficients Cj can 
be obtained systematically for any m. Fig. |2l illustrates the 
convergence of the total energy with respect to the number 
of grid points using 5-point (m — 2), 13-point (m — 6) 
FD and PiB representations. The lower-order finite difference 
scheme, m = 2, is poor in terms of accuracy, and converges 
slowly as the grid size increases. The high-order FD scheme, 
TO = 6, does improve the accuracy by two orders of magni- 
tude, but is still less accurate than PiB for Nr^ = 32, 48 and 
64. 

We check the accuracy and efficiency of our line- 
minimization scheme by comparing with TPA's approach as 
well as the numerically exact Brent's line search algorithm,^*' 
In this case, we use A^band = 5 and iVupdate = 1 as recom- 
mended in TPA's original algorithm.^' Fig. |3 plots in 
the first 25 band iterations calculated from three schemes re- 
spectively. The values of 6',ni„ calculated from both TPA's and 
our method agree very well with exact values and the relative 
errors are always smaller than 1%. But in terms of computa- 
tional effort, our new scheme is much more efficient as argued 
in the previous section. 

To find the optimal A'^band and iVupdatc, we do the calcu- 
lations with different values of iVband and iVupdatc, and the 
results are shown in Fig. |3 With fixed A^updatc = 1, it is 
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FIG. 3: Comparison of ^min calculated by three different schemes. 
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FIG. 4: Errors in the total energy as a function of CPU time during 
the DMCG calculation for different A^band with A^'updatc = 1 (a) 
and for different TVupdatc with A'^band ~ 20 (b) with respect to the 
exact total energy that is calculated using a finer grid (Nx = 80) and 
tighter convergence criterion (e = 10~^) , 



seen that a relatively larger A'^band is more efficient than small 
iVband- Fixing iVband = 20, A^updatc = 20 gives the best per- 
formance. The combination of large A'^band and A/updatc re- 
duces the computational effort by almost one order of magni- 
tude. Though the actual values of optimal A/band and iVupdatc 
may vary for different systems, the basic idea demonstrated in 
this test calculation is believed to be of general significance. 

We have implemented both the two-level (h and 2h) and 
three-level and 2h) SMG schemes. Instead of using a 

sophisticated interpolation as in Ref. |33l we use a simpler 
Lagrange polynomial interpolation method. Comparison with 
the single-level calculation shows a quite obvious improve- 
ment in computational efficiency. To check the effect of the 
interpolation accuracy, in Fig. |5lwe plot the relative compu- 



tational effort in one KS calculation as a function of the order 
of the Lagrange interpolation formula in both two- and three- 
level SMG calculations. We see that a high-order interpola- 
tion scheme is useful to improve the performance. 




Two-Level 
Three- Level 



2 4 6 8 

Order of Lagrange Interpolation 

FIG. 5: Relative computational efforts for one KS calculation as 
a function of the order of interpolation in both two- and three-level 
SMG calculations. The CPU time in the single-level calculation is 
taken as the unit. 



IV. SUMMARY 

In this paper, we presented an efficient algorithm for the 
KS-SDFT simulation of large quantum dot systems. The main 
elements of the algorithm are: (1) Wave functions are repre- 
sented in real space, and the kinetic energy operator is applied 
to wave functions by fast sine transform. (2) The Hartree po- 
tential is calculated by Martyna-Tuckerman's Fourier convo- 
lution method. (3) For the solution of the KS equation, we in- 
troduced several important modifications to Teter et al.'s band- 
by-band conjugate-gradient method. A more efficient approx- 
imate line-minimization scheme was developed; it was found 
that large band iteration number and a delayed update of the 
KS Hamiltonian inside the band iterations increase the effi- 
ciency by one order of magnitude. (4) A simplified multi-grid 
technique was introduced to accelerate the convergence. The 
new algorithm has been used to study spin and conductance 
peak-spacing distributions in a 2D chaotic QD system with 
electron number N up to 200, from which new physical phe- 
nomena were revealedii^ 
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